%Script to create a mask to omit US states that do not border any of the Great Lakes (for visualization purposes)

clear

addpath(genpath('/home/mmfire/matlab_functions/FUNCTIONS'))
addpath(genpath('/home/mmfire/matlab_functions/SHAPEFILES'))

CurrentFolder=pwd;

dirp=strcat('/Volumes/TOSHIBA/climgroup_research_macbookair/HWA/manuscripts/C__HWAmicroclimate/companion_dataset/05_prism/');%prism [AN81d dataset;12Z-12Z PRISM day]

%Contents of *bil.hdr file.
BYTEORDER='I';
LAYOUT='BIL';
NROWS=621;
NCOLS=1405;
NBANDS=1;
NBITS=32;
BANDROWBYTES=5620;
TOTALROWBYTES=5620;
PIXELTYPE='FLOAT';
ULXMAP=-125;
ULYMAP=49.9166666666664;
XDIM=0.0416666666667;
YDIM=0.0416666666667;
NODATA=-9999;

%Construct latitude and longitude arrays
%From : http://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/bil-bip-and-bsq-raster-files.htm
%ulxmap—The x-axis map coordinate of the center of the upper left pixel. If you specify this parameter, set ulymap, too; otherwise, a default value is used.
%ulymap—The y-axis map coordinate of the center of the upper left pixel. If this parameter is specified, ulxmap must also be set; otherwise, a default value is used.
temLATS=nan(NROWS,NCOLS);temLONS=nan(NROWS,NCOLS);
for j=NROWS:-1:1
    temLATS(j,:)=ULYMAP+((j-NROWS)*YDIM);
end
for i=1:1:NCOLS
    temLONS(:,i)=ULXMAP+((i-1)*XDIM);
end
%Subset full PRISM CONUS domain: this is approximately equal to B.J.'s
%subdomain
ss_w=-98.6;
ss_e=-71.3;
ss_s=36.9;
ss_n=49.6;
I=(temLATS>=ss_s)&(temLATS<=ss_n)&(temLONS>=ss_w)&(temLONS<=ss_e);
[jll,ill] = find(I>0,1,'first');[jur,iur] = find(I>0,1,'last');clear I

LATS=temLATS(jll:jur,ill:iur);LONS=temLONS(jll:jur,ill:iur);

%Create multistate mask: 
%Mask out values outside of: MN,WI,MI,IL,IN,OH,PA,NY
mask=LATS*0;
HWApoly = shaperead('usastatehi', 'UseGeoCoords', true,...
  'Selector',...
  {@(name) any(strcmp(name,{'Minnesota'}) | strcmp(name,{'Wisconsin'})...
  | strcmp(name,{'Michigan'}) | strcmp(name,{'Illinois'})...
  | strcmp(name,{'Indiana'}) | strcmp(name,{'Ohio'})...
  | strcmp(name,{'Pennsylvania'}) | strcmp(name,{'New York'})), 'Name'});
for j=1:size(LATS,1)
    for i=1:size(LONS,2)
        disp(strcat(num2str(i),',',num2str(j)))
        c=0;
        for s=1:size(HWApoly,1)
            if(c==0)
                SLAT=HWApoly(s).Lat;SLON=HWApoly(s).Lon;
                I=inpolygon(LONS(j,i),LATS(j,i),SLON,SLAT);
                if(I==1);mask(j,i) = I;c=1;end
            end
        end
    end
end
save('hwastatemask.mat','mask')
